﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides the numerical integration for integrals by using the Romberg 
    /// algorithm.
    /// </summary>
    [Serializable]
    public class RombergIntegrator : AbstractCompositeIntegrator
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="RombergIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        public RombergIntegrator(double a, double b, IRealFunction function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="RombergIntegrator"/> class.
        /// </summary>
        /// <param name="a">A.</param>
        /// <param name="b">The b.</param>
        /// <param name="function">The function.</param>
        public RombergIntegrator(double a, double b, SimplePolynomial function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="RombergIntegrator"/> class.
        /// </summary>
        /// <param name="a">A.</param>
        /// <param name="b">The b.</param>
        /// <param name="function">The function.</param>
        public RombergIntegrator(double a, double b, Polynomial function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="RombergIntegrator"/> class.
        /// </summary>
        /// <param name="a">A.</param>
        /// <param name="b">The b.</param>
        /// <param name="function">The function.</param>
        public RombergIntegrator(double a, double b, SimpleRationalPolynomial function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="RombergIntegrator"/> class.
        /// </summary>
        /// <param name="a">A.</param>
        /// <param name="b">The b.</param>
        /// <param name="function">The function.</param>
        public RombergIntegrator(double a, double b, HardRealFunction function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="RombergIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The integral for the numerical integration.</param>
        public RombergIntegrator(FinitelyDoubleIntegral integral)
            : base(integral.A, integral.B, integral.Function)
        {
        }

        /// <summary>
        /// Integrates the integral in a numerical way. After running the numerical integration
        /// it should be checked if the result has a precision error.
        /// </summary>
        /// <param name="iterations">The maximum numbers of iterations.</param>
        /// <param name="precision">The minimum precision of the numerical integration.</param>
        /// <returns>The area of the specified integral.</returns>
        public override double Integrate(int iterations, double precision)
        {
            bool swapBounds = false;

            if (this.B < this.A)
            {
                double tempuri = this.A;

                this.A = this.B;
                this.B = tempuri;
                swapBounds = true;
            }

            int i = 1;
            double[,] workingData = new double[iterations,iterations];

            workingData[0, 0] = 0.5 * (this.B - this.A) * (this.Function.SolveAt(this.A)
                                                           + this.Function.SolveAt(this.B));

            while (i < iterations)
            {
                double tempuri = 0;
                double h = (this.B - this.A) / Math.Pow(2.0, i);

                for (int j = 0; j <= Math.Pow(2, i - 1); j++)
                {
                    tempuri += (this.Function.SolveAt(this.A + (2 * j) * h));
                }

                workingData[i, 0] = 0.5 * workingData[i - 1, 0] + h * tempuri;

                for (int j = 0; j <= i; j++)
                {
                    workingData[i, j + 1] = workingData[i, j] +
                                            (workingData[i, j] - workingData[i - 1, j]) / (4 * (j + 1) - 1);
                }

                this.relativeError = Math.Abs(workingData[i, i - 1] - workingData[i, i]);

                if (this.relativeError < precision)
                {
                    this.neededIterations = i;
                    this.precisionError = false;

                    if (swapBounds)
                    {
                        double tempuri1 = this.A;

                        this.A = this.B;
                        this.B = tempuri1;

                        return -workingData[i, i];
                    }

                    return workingData[i, i];
                }

                i++;
            }

            this.neededIterations = i;
            this.precisionError = true;

            if (swapBounds)
            {
                double tempuri1 = this.A;

                this.A = this.B;
                this.B = tempuri1;

                return -workingData[i, i];
            }

            return workingData[i, i];
        }
    }
}